#!/usr/bin/env python
# -*- coding: utf-8 -*-

""" By Martin Senande-Rivera
    For Towards and atmosphere more favourable to firestorm development in Europe """

import os
import glob, sys
import numpy as np
import xarray as xr
import pandas as pd
from scipy import stats

GCM=sys.argv[1]  # Enter GCM
RCM=sys.argv[2]  # Enter RCM
rcp=sys.argv[3]  # Enter RCP scenario

path_outs='./'+RCM+'/'+GCM+'/'

# Number of days with Kindex above the threshold
Ki = xr.open_dataset(path_outs+'Ki_ndays_'+rcp+'.nc')['__xarray_dataarray_variable__']
Ki = Ki.sel(time=slice('1970-01-01','2098-12-31'))

# Linear trend function
def linear_trend(x):
    pf = np.polyfit(np.arange(1970.,2099.), x, 1)
    return xr.DataArray(pf[0])
# p-values function
def pvalue(x):
    r, p = stats.pearsonr(np.arange(1970.,2099.), x)
    return xr.DataArray(p)

# Linear trend of the number of days with Kindex above the threshold
trend = xr.apply_ufunc(linear_trend,  # function to apply
                         Ki,
                         input_core_dims=(('time',),),
                         vectorize=True) 
trend.to_netcdf(path_outs+'Ki_ndays_trend_'+rcp+'.nc')

# p-values of the number of days with Kindex above the threshold
pvalues = xr.apply_ufunc(pvalue,  # function to apply
                         Ki,
                         input_core_dims=(('time',),),
                         vectorize=True) 
pvalues.to_netcdf(path_outs+'Ki_ndays_pvalues_'+rcp+'.nc')

